home *** CD-ROM | disk | FTP | other *** search
/ BCI NET 2 / BCI NET 2.iso / archives / programming / source / graphicgems4.lha / GemsIV / ptpoly_haines / ptinpoly.c < prev    next >
Encoding:
C/C++ Source or Header  |  1995-02-06  |  52.3 KB  |  2,046 lines
  1. /* ptinpoly.c - point in polygon inside/outside code.
  2.  
  3.    by Eric Haines, 3D/Eye Inc, erich@eye.com
  4.  
  5.    This code contains the following algorithms:
  6.     crossings - count the crossing made by a ray from the test point
  7.     angle summation - sum the angle formed by point and vertex pairs
  8.     weiler angle summation - sum the angles using quad movements
  9.     half-plane testing - test triangle fan using half-space planes
  10.     barycentric coordinates - test triangle fan w/barycentric coords
  11.     spackman barycentric - preprocessed barycentric coordinates
  12.     trapezoid testing - bin sorting algorithm
  13.     grid testing - grid imposed on polygon
  14.     exterior test - for convex polygons, check exterior of polygon
  15.     inclusion test - for convex polygons, use binary search for edge.
  16. */
  17.  
  18. #include <stdio.h>
  19. #include <stdlib.h>
  20. #include <math.h>
  21. #include "ptinpoly.h"
  22.  
  23. #define X    0
  24. #define Y    1
  25.  
  26. #ifndef TRUE
  27. #define TRUE    1
  28. #define FALSE    0
  29. #endif
  30.  
  31. #ifndef HUGE
  32. #define HUGE    1.79769313486232e+308
  33. #endif
  34.  
  35. /* test if a & b are within epsilon.  Favors cases where a < b */
  36. #define Near(a,b,eps)    ( ((b)-(eps)<(a)) && ((a)-(eps)<(b)) )
  37.  
  38. #define MALLOC_CHECK( a )    if ( !(a) ) {                   \
  39.                     fprintf( stderr, "out of memory\n" ) ; \
  40.                     exit(1) ;                   \
  41.                 }
  42.  
  43.  
  44. /* ======= Crossings algorithm ============================================ */
  45.  
  46. /* Shoot a test ray along +X axis.  The strategy, from MacMartin, is to
  47.  * compare vertex Y values to the testing point's Y and quickly discard
  48.  * edges which are entirely to one side of the test ray.
  49.  *
  50.  * Input 2D polygon _pgon_ with _numverts_ number of vertices and test point
  51.  * _point_, returns 1 if inside, 0 if outside.    WINDING and CONVEX can be
  52.  * defined for this test.
  53.  */
  54. int CrossingsTest( pgon, numverts, point )
  55. double    pgon[][2] ;
  56. int    numverts ;
  57. double    point[2] ;
  58. {
  59. #ifdef    WINDING
  60. register int    crossings ;
  61. #endif
  62. register int    j, yflag0, yflag1, inside_flag, xflag0 ;
  63. register double ty, tx, *vtx0, *vtx1 ;
  64. #ifdef    CONVEX
  65. register int    line_flag ;
  66. #endif
  67.  
  68.     tx = point[X] ;
  69.     ty = point[Y] ;
  70.  
  71.     vtx0 = pgon[numverts-1] ;
  72.     /* get test bit for above/below X axis */
  73.     yflag0 = ( vtx0[Y] >= ty ) ;
  74.     vtx1 = pgon[0] ;
  75.  
  76. #ifdef    WINDING
  77.     crossings = 0 ;
  78. #else
  79.     inside_flag = 0 ;
  80. #endif
  81. #ifdef    CONVEX
  82.     line_flag = 0 ;
  83. #endif
  84.     for ( j = numverts+1 ; --j ; ) {
  85.  
  86.     yflag1 = ( vtx1[Y] >= ty ) ;
  87.     /* check if endpoints straddle (are on opposite sides) of X axis
  88.      * (i.e. the Y's differ); if so, +X ray could intersect this edge.
  89.      */
  90.     if ( yflag0 != yflag1 ) {
  91.         xflag0 = ( vtx0[X] >= tx ) ;
  92.         /* check if endpoints are on same side of the Y axis (i.e. X's
  93.          * are the same); if so, it's easy to test if edge hits or misses.
  94.          */
  95.         if ( xflag0 == ( vtx1[X] >= tx ) ) {
  96.  
  97.         /* if edge's X values both right of the point, must hit */
  98. #ifdef    WINDING
  99.         if ( xflag0 ) crossings += ( yflag0 ? -1 : 1 ) ;
  100. #else
  101.         if ( xflag0 ) inside_flag = !inside_flag ;
  102. #endif
  103.         } else {
  104.         /* compute intersection of pgon segment with +X ray, note
  105.          * if >= point's X; if so, the ray hits it.
  106.          */
  107.         if ( (vtx1[X] - (vtx1[Y]-ty)*
  108.              ( vtx0[X]-vtx1[X])/(vtx0[Y]-vtx1[Y])) >= tx ) {
  109. #ifdef    WINDING
  110.             crossings += ( yflag0 ? -1 : 1 ) ;
  111. #else
  112.             inside_flag = !inside_flag ;
  113. #endif
  114.         }
  115.         }
  116. #ifdef    CONVEX
  117.         /* if this is second edge hit, then done testing */
  118.         if ( line_flag ) goto Exit ;
  119.  
  120.         /* note that one edge has been hit by the ray's line */
  121.         line_flag = TRUE ;
  122. #endif
  123.     }
  124.  
  125.     /* move to next pair of vertices, retaining info as possible */
  126.     yflag0 = yflag1 ;
  127.     vtx0 = vtx1 ;
  128.     vtx1 += 2 ;
  129.     }
  130. #ifdef    CONVEX
  131.     Exit: ;
  132. #endif
  133. #ifdef    WINDING
  134.     /* test if crossings is not zero */
  135.     inside_flag = (crossings != 0) ;
  136. #endif
  137.  
  138.     return( inside_flag ) ;
  139. }
  140.  
  141. /* ======= Angle summation algorithm ======================================= */
  142.  
  143. /* Sum angles made by (vtxN to test point to vtxN+1), check if in proper
  144.  * range to be inside.    VERY SLOW, included for tutorial reasons (i.e.
  145.  * to show why one should never use this algorithm).
  146.  *
  147.  * Input 2D polygon _pgon_ with _numverts_ number of vertices and test point
  148.  * _point_, returns 1 if inside, 0 if outside.
  149.  */
  150. int AngleTest( pgon, numverts, point )
  151. double    pgon[][2] ;
  152. int    numverts ;
  153. double    point[2] ;
  154. {
  155. register double *vtx0, *vtx1, angle, len, vec0[2], vec1[2], vec_dot ;
  156. register int    j ;
  157. int    inside_flag ;
  158.  
  159.     /* sum the angles and see if answer mod 2*PI > PI */
  160.     vtx0 = pgon[numverts-1] ;
  161.     vec0[X] = vtx0[X] - point[X] ;
  162.     vec0[Y] = vtx0[Y] - point[Y] ;
  163.     len = sqrt( vec0[X] * vec0[X] + vec0[Y] * vec0[Y] ) ;
  164.     if ( len <= 0.0 ) {
  165.     /* point and vertex coincide */
  166.     return( 1 ) ;
  167.     }
  168.     vec0[X] /= len ;
  169.     vec0[Y] /= len ;
  170.  
  171.     angle = 0.0 ;
  172.     for ( j = 0 ; j < numverts ; j++ ) {
  173.     vtx1 = pgon[j] ;
  174.     vec1[X] = vtx1[X] - point[X] ;
  175.     vec1[Y] = vtx1[Y] - point[Y] ;
  176.     len = sqrt( vec1[X] * vec1[X] + vec1[Y] * vec1[Y] ) ;
  177.     if ( len <= 0.0 ) {
  178.         /* point and vertex coincide */
  179.         return( 1 ) ;
  180.     }
  181.     vec1[X] /= len ;
  182.     vec1[Y] /= len ;
  183.  
  184.     /* check if vec1 is to "left" or "right" of vec0 */
  185.     vec_dot = vec0[X] * vec1[X] + vec0[Y] * vec1[Y] ;
  186.     if ( vec_dot < -1.0 ) {
  187.         /* point is on edge, so always add 180 degrees.  always
  188.          * adding is not necessarily the right answer for
  189.          * concave polygons and subtractive triangles.
  190.          */
  191.         angle += M_PI ;
  192.     } else if ( vec_dot < 1.0 ) {
  193.         if ( vec0[X] * vec1[Y] - vec1[X] * vec0[Y] >= 0.0 ) {
  194.         /* add angle due to dot product of vectors */
  195.         angle += acos( vec_dot ) ;
  196.         } else {
  197.         /* subtract angle due to dot product of vectors */
  198.         angle -= acos( vec_dot ) ;
  199.         }
  200.     } /* if vec_dot >= 1.0, angle does not change */
  201.  
  202.     /* get to next point */
  203.     vtx0 = vtx1 ;
  204.     vec0[X] = vec1[X] ;
  205.     vec0[Y] = vec1[Y] ;
  206.     }
  207.     /* test if between PI and 3*PI, 5*PI and 7*PI, etc */
  208.     inside_flag = fmod( fabs(angle) + M_PI, 4.0*M_PI ) > 2.0*M_PI ;
  209.  
  210.     return( inside_flag ) ;
  211. }
  212.  
  213. /* ======= Weiler algorithm ============================================ */
  214.  
  215. /* Track quadrant movements using Weiler's algorithm (elsewhere in Graphics
  216.  * Gems IV).  Algorithm has been optimized for testing purposes, but the
  217.  * crossings algorithm is still faster.     Included to provide timings.
  218.  *
  219.  * Input 2D polygon _pgon_ with _numverts_ number of vertices and test point
  220.  * _point_, returns 1 if inside, 0 if outside.    WINDING can be defined for
  221.  * this test.
  222.  */
  223.  
  224. #define QUADRANT( vtx, x, y )    \
  225.     (((vtx)[X]>(x)) ? ( ((vtx)[Y]>(y)) ? 0:3 ) : ( ((vtx)[Y]>(y)) ? 1:2 ))
  226.  
  227. #define X_INTERCEPT( v0, v1, y )    \
  228.     ( (((v1)[X]-(v0)[X])/((v1)[Y]-(v0)[Y])) * ((y)-(v0)[Y]) + (v0)[X] )
  229.  
  230. int WeilerTest( pgon, numverts, point )
  231. double    pgon[][2] ;
  232. int    numverts ;
  233. double    point[2] ;
  234. {
  235. register int    angle, qd, next_qd, delta, j ;
  236. register double ty, tx, *vtx0, *vtx1 ;
  237. int    inside_flag ;
  238.  
  239.     tx = point[X] ;
  240.     ty = point[Y] ;
  241.  
  242.     vtx0 = pgon[numverts-1] ;
  243.     qd = QUADRANT( vtx0, tx, ty ) ;
  244.     angle = 0 ;
  245.  
  246.     vtx1 = pgon[0] ;
  247.  
  248.     for ( j = numverts+1 ; --j ; ) {
  249.     /* calculate quadrant and delta from last quadrant */
  250.     next_qd = QUADRANT( vtx1, tx, ty ) ;
  251.     delta = next_qd - qd ;
  252.  
  253.     /* adjust delta and add it to total angle sum */
  254.     switch ( delta ) {
  255.         case 0:    /* do nothing */
  256.         break ;
  257.         case -1:
  258.         case 3:
  259.         angle-- ;
  260.         qd = next_qd ;
  261.         break ;
  262.  
  263.         case 1:
  264.         case -3:
  265.         angle++ ;
  266.         qd = next_qd ;
  267.         break ;
  268.  
  269.         case 2:
  270.         case -2:
  271.         if (X_INTERCEPT( vtx0, vtx1, ty ) > tx ) {
  272.             angle -= delta ;
  273.         } else {
  274.             angle += delta ;
  275.         }
  276.         qd = next_qd ;
  277.         break ;
  278.     }
  279.  
  280.     /* increment for next step */
  281.     vtx0 = vtx1 ;
  282.     vtx1 += 2 ;
  283.     }
  284.  
  285. #ifdef    WINDING
  286.     /* simple windings test:  if angle != 0, then point is inside */
  287.     inside_flag = ( angle != 0 ) ;
  288. #else
  289.     /* Jordan test:  if angle is +-4, 12, 20, etc then point is inside */
  290.     inside_flag = ( (angle/4) & 0x1 ) ;
  291. #endif
  292.  
  293.     return (inside_flag) ;
  294. }
  295.  
  296. #undef    QUADRANT
  297. #undef    Y_INTERCEPT
  298.  
  299. /* ======= Triangle half-plane algorithm ================================== */
  300.  
  301. /* Split the polygon into a fan of triangles and for each triangle test if
  302.  * the point is inside of the three half-planes formed by the triangle's edges.
  303.  *
  304.  * Call setup with 2D polygon _pgon_ with _numverts_ number of vertices,
  305.  * which returns a pointer to a plane set array.
  306.  * Call testing procedure with a pointer to this array, _numverts_, and
  307.  * test point _point_, returns 1 if inside, 0 if outside.
  308.  * Call cleanup with pointer to plane set array to free space.
  309.  *
  310.  * SORT and CONVEX can be defined for this test.
  311.  */
  312.  
  313. /* split polygons along set of x axes - call preprocess once */
  314. pPlaneSet    PlaneSetup( pgon, numverts )
  315. double    pgon[][2] ;
  316. int    numverts ;
  317. {
  318. int    i, p1, p2 ;
  319. double    tx, ty, vx0, vy0 ;
  320. pPlaneSet    pps, pps_return ;
  321. #ifdef    SORT
  322. double    len[3], len_temp ;
  323. int    j ;
  324. PlaneSet    ps_temp ;
  325. #ifdef    CONVEX
  326. pPlaneSet    pps_new ;
  327. pSizePlanePair p_size_pair ;
  328. #endif
  329. #endif
  330.  
  331.     pps = pps_return =
  332.         (pPlaneSet)malloc( 3 * (numverts-2) * sizeof( PlaneSet )) ;
  333.     MALLOC_CHECK( pps ) ;
  334. #ifdef    CONVEX
  335. #ifdef    SORT
  336.     p_size_pair =
  337.     (pSizePlanePair)malloc( (numverts-2) * sizeof( SizePlanePair ) ) ;
  338.     MALLOC_CHECK( p_size_pair ) ;
  339. #endif
  340. #endif
  341.  
  342.     vx0 = pgon[0][X] ;
  343.     vy0 = pgon[0][Y] ;
  344.  
  345.     for ( p1 = 1, p2 = 2 ; p2 < numverts ; p1++, p2++ ) {
  346.     pps->vx = vy0 - pgon[p1][Y] ;
  347.     pps->vy = pgon[p1][X] - vx0 ;
  348.     pps->c = pps->vx * vx0 + pps->vy * vy0 ;
  349. #ifdef    SORT
  350.     len[0] = pps->vx * pps->vx + pps->vy * pps->vy ;
  351. #ifdef    CONVEX
  352. #ifdef    HYBRID
  353.     pps->ext_flag = ( p1 == 1 ) ;
  354. #endif
  355.     /* Sort triangles by areas, so compute (twice) the area here */
  356.     p_size_pair[p1-1].pps = pps ;
  357.     p_size_pair[p1-1].size =
  358.             ( pgon[0][X] * pgon[p1][Y] ) +
  359.             ( pgon[p1][X] * pgon[p2][Y] ) +
  360.             ( pgon[p2][X] * pgon[0][Y] ) -
  361.             ( pgon[p1][X] * pgon[0][Y] ) -
  362.             ( pgon[p2][X] * pgon[p1][Y] ) -
  363.             ( pgon[0][X] * pgon[p2][Y] ) ;
  364. #endif
  365. #endif
  366.     pps++ ;
  367.     pps->vx = pgon[p1][Y] - pgon[p2][Y] ;
  368.     pps->vy = pgon[p2][X] - pgon[p1][X] ;
  369.     pps->c = pps->vx * pgon[p1][X] + pps->vy * pgon[p1][Y] ;
  370. #ifdef    SORT
  371.     len[1] = pps->vx * pps->vx + pps->vy * pps->vy ;
  372. #endif
  373. #ifdef    CONVEX
  374. #ifdef    HYBRID
  375.     pps->ext_flag = TRUE ;
  376. #endif
  377. #endif
  378.     pps++ ;
  379.     pps->vx = pgon[p2][Y] - vy0 ;
  380.     pps->vy = vx0 - pgon[p2][X] ;
  381.     pps->c = pps->vx * pgon[p2][X] + pps->vy * pgon[p2][Y] ;
  382. #ifdef    SORT
  383.     len[2] = pps->vx * pps->vx + pps->vy * pps->vy ;
  384. #endif
  385. #ifdef    CONVEX
  386. #ifdef    HYBRID
  387.     pps->ext_flag = ( p2 == numverts-1 ) ;
  388. #endif
  389. #endif
  390.  
  391.     /* find an average point which must be inside of the triangle */
  392.     tx = ( vx0 + pgon[p1][X] + pgon[p2][X] ) / 3.0 ;
  393.     ty = ( vy0 + pgon[p1][Y] + pgon[p2][Y] ) / 3.0 ;
  394.  
  395.     /* check sense and reverse if test point is not thought to be inside
  396.      * first triangle
  397.      */
  398.     if ( pps->vx * tx + pps->vy * ty >= pps->c ) {
  399.         /* back up to start of plane set */
  400.         pps -= 2 ;
  401.         /* point is thought to be outside, so reverse sense of edge
  402.          * normals so that it is correctly considered inside.
  403.          */
  404.         for ( i = 0 ; i < 3 ; i++ ) {
  405.         pps->vx = -pps->vx ;
  406.         pps->vy = -pps->vy ;
  407.         pps->c    = -pps->c ;
  408.         pps++ ;
  409.         }
  410.     } else {
  411.         pps++ ;
  412.     }
  413.  
  414. #ifdef    SORT
  415.     /* sort the planes based on the edge lengths */
  416.     pps -= 3 ;
  417.     for ( i = 0 ; i < 2 ; i++ ) {
  418.         for ( j = i+1 ; j < 3 ; j++ ) {
  419.         if ( len[i] < len[j] ) {
  420.             ps_temp = pps[i] ;
  421.             pps[i] = pps[j] ;
  422.             pps[j] = ps_temp ;
  423.             len_temp = len[i] ;
  424.             len[i] = len[j] ;
  425.             len[j] = len_temp ;
  426.         }
  427.         }
  428.     }
  429.     pps += 3 ;
  430. #endif
  431.     }
  432.  
  433. #ifdef    CONVEX
  434. #ifdef    SORT
  435.     /* sort the triangles based on their areas */
  436.     qsort( p_size_pair, numverts-2,
  437.         sizeof( SizePlanePair ), CompareSizePlanePairs ) ;
  438.  
  439.     /* make the plane sets match the sorted order */
  440.     for ( i = 0, pps = pps_return
  441.     ; i < numverts-2
  442.     ; i++ ) {
  443.  
  444.     pps_new = p_size_pair[i].pps ;
  445.     for ( j = 0 ; j < 3 ; j++, pps++, pps_new++ ) {
  446.         ps_temp = *pps ;
  447.         *pps = *pps_new ;
  448.         *pps_new = ps_temp ;
  449.     }
  450.     }
  451.     free( p_size_pair ) ;
  452. #endif
  453. #endif
  454.  
  455.     return( pps_return ) ;
  456. }
  457.  
  458. #ifdef    CONVEX
  459. #ifdef    SORT
  460. int CompareSizePlanePairs( p_sp0, p_sp1 )
  461. pSizePlanePair    p_sp0, p_sp1 ;
  462. {
  463.     if ( p_sp0->size == p_sp1->size ) {
  464.     return( 0 ) ;
  465.     } else {
  466.     return( p_sp0->size > p_sp1->size ? -1 : 1 ) ;
  467.     }
  468. }
  469. #endif
  470. #endif
  471.  
  472.  
  473. /* check point for inside of three "planes" formed by triangle edges */
  474. int PlaneTest( p_plane_set, numverts, point )
  475. pPlaneSet    p_plane_set ;
  476. int    numverts ;
  477. double    point[2] ;
  478. {
  479. register pPlaneSet    ps ;
  480. register int    p2 ;
  481. #ifndef CONVEX
  482. register int    inside_flag ;
  483. #endif
  484. register double tx, ty ;
  485.  
  486.     tx = point[X] ;
  487.     ty = point[Y] ;
  488.  
  489. #ifndef CONVEX
  490.     inside_flag = 0 ;
  491. #endif
  492.  
  493.     for ( ps = p_plane_set, p2 = numverts-1 ; --p2 ; ) {
  494.  
  495.     if ( ps->vx * tx + ps->vy * ty < ps->c ) {
  496.         ps++ ;
  497.         if ( ps->vx * tx + ps->vy * ty < ps->c ) {
  498.         ps++ ;
  499.         /* note: we make the third edge have a slightly different
  500.          * equality condition, since this third edge is in fact
  501.          * the next triangle's first edge.  Not fool-proof, but
  502.          * it doesn't hurt (better would be to keep track of the
  503.          * triangle's area sign so we would know which kind of
  504.          * triangle this is).  Note that edge sorting nullifies
  505.          * this special inequality, too.
  506.          */
  507.         if ( ps->vx * tx + ps->vy * ty <= ps->c ) {
  508.             /* point is inside polygon */
  509. #ifdef CONVEX
  510.             return( 1 ) ;
  511. #else
  512.             inside_flag = !inside_flag ;
  513. #endif
  514.         }
  515. #ifdef    CONVEX
  516. #ifdef    HYBRID
  517.         /* check if outside exterior edge */
  518.         else if ( ps->ext_flag ) return( 0 ) ;
  519. #endif
  520. #endif
  521.         ps++ ;
  522.         } else {
  523. #ifdef    CONVEX
  524. #ifdef    HYBRID
  525.         /* check if outside exterior edge */
  526.         if ( ps->ext_flag ) return( 0 ) ;
  527. #endif
  528. #endif
  529.         /* get past last two plane tests */
  530.         ps += 2 ;
  531.         }
  532.     } else {
  533. #ifdef    CONVEX
  534. #ifdef    HYBRID
  535.         /* check if outside exterior edge */
  536.         if ( ps->ext_flag ) return( 0 ) ;
  537. #endif
  538. #endif
  539.         /* get past all three plane tests */
  540.         ps += 3 ;
  541.     }
  542.     }
  543.  
  544. #ifdef CONVEX
  545.     /* for convex, if we make it to here, all triangles were missed */
  546.     return( 0 ) ;
  547. #else
  548.     return( inside_flag ) ;
  549. #endif
  550. }
  551.  
  552. void PlaneCleanup( p_plane_set )
  553. pPlaneSet    p_plane_set ;
  554. {
  555.     free( p_plane_set ) ;
  556. }
  557.  
  558. /* ======= Barycentric algorithm ========================================== */
  559.  
  560. /* Split the polygon into a fan of triangles and for each triangle test if
  561.  * the point has barycentric coordinates which are inside the triangle.
  562.  * Similar to Badouel's code in Graphics Gems, with a little more efficient
  563.  * coding.
  564.  *
  565.  * Input 2D polygon _pgon_ with _numverts_ number of vertices and test point
  566.  * _point_, returns 1 if inside, 0 if outside.
  567.  */
  568. int BarycentricTest( pgon, numverts, point )
  569. double    pgon[][2] ;
  570. int    numverts ;
  571. double    point[2] ;
  572. {
  573. register double *pg1, *pg2, *pgend ;
  574. register double tx, ty, u0, u1, u2, v0, v1, vx0, vy0, alpha, beta, denom ;
  575. int    inside_flag ;
  576.  
  577.     tx = point[X] ;
  578.     ty = point[Y] ;
  579.     vx0 = pgon[0][X] ;
  580.     vy0 = pgon[0][Y] ;
  581.     u0 = tx - vx0 ;
  582.     v0 = ty - vy0 ;
  583.  
  584.     inside_flag = 0 ;
  585.     pgend = pgon[numverts-1] ;
  586.     for ( pg1 = pgon[1], pg2 = pgon[2] ; pg1 != pgend ; pg1+=2, pg2+=2 ) {
  587.  
  588.     u1 = pg1[X] - vx0 ;
  589.     if ( u1 == 0.0 ) {
  590.  
  591.         /* 0 and 1 vertices have same X value */
  592.  
  593.         /* zero area test - can be removed for convex testing */
  594.         u2 = pg2[X] - vx0 ;
  595.         if ( ( u2 == 0.0 ) ||
  596.  
  597.         /* compute beta and check bounds */
  598.         /* we use "<= 0.0" so that points on the shared interior
  599.          * edge will (generally) be inside only one polygon.
  600.          */
  601.          ( ( beta = u0 / u2 ) <= 0.0 ) ||
  602.          ( beta > 1.0 ) ||
  603.  
  604.          /* zero area test - remove for convex testing */
  605.          ( ( v1 = pg1[Y] - vy0 ) == 0.0 ) ||
  606.  
  607.          /* compute alpha and check bounds */
  608.          ( ( alpha = ( v0 - beta *
  609.             ( pg2[Y] - vy0 ) ) / v1 ) < 0.0 ) ) {
  610.  
  611.         /* whew! missed! */
  612.         goto NextTri ;
  613.         }
  614.  
  615.     } else {
  616.         /* 0 and 1 vertices have different X value */
  617.  
  618.         /* compute denom and check for zero area triangle - check
  619.          * is not needed for convex polygon testing
  620.          */
  621.         u2 = pg2[X] - vx0 ;
  622.         v1 = pg1[Y] - vy0 ;
  623.         denom = ( pg2[Y] - vy0 ) * u1 - u2 * v1 ;
  624.         if ( ( denom == 0.0 ) ||
  625.  
  626.         /* compute beta and check bounds */
  627.         /* we use "<= 0.0" so that points on the shared interior
  628.          * edge will (generally) be inside only one polygon.
  629.          */
  630.          ( ( beta = ( v0 * u1 - u0 * v1 ) / denom ) <= 0.0 ) ||
  631.          ( beta > 1.0 ) ||
  632.  
  633.          /* compute alpha & check bounds */
  634.          ( ( alpha = ( u0 - beta * u2 ) / u1 ) < 0.0 ) ) {
  635.  
  636.         /* whew! missed! */
  637.         goto NextTri ;
  638.         }
  639.     }
  640.  
  641.     /* check gamma */
  642.     if ( alpha + beta <= 1.0 ) {
  643.         /* survived */
  644.         inside_flag = !inside_flag ;
  645.     }
  646.  
  647.     NextTri: ;
  648.     }
  649.     return( inside_flag ) ;
  650. }
  651.  
  652. /* ======= Barycentric precompute (Spackman) algorithm ===================== */
  653.  
  654. /* Split the polygon into a fan of triangles and for each triangle test if
  655.  * the point has barycentric coordinates which are inside the triangle.
  656.  * Use Spackman's normalization method to precompute various parameters.
  657.  *
  658.  * Call setup with 2D polygon _pgon_ with _numverts_ number of vertices,
  659.  * which returns a pointer to the array of the parameters records and the
  660.  * number of parameter records created.
  661.  * Call testing procedure with the first vertex in the polygon _pgon[0]_,
  662.  * a pointer to this array, the number of parameter records, and test point
  663.  * _point_, returns 1 if inside, 0 if outside.
  664.  * Call cleanup with pointer to parameter record array to free space.
  665.  *
  666.  * SORT can be defined for this test.
  667.  * (CONVEX could be added: see PlaneSetup and PlaneTest for method)
  668.  */
  669. pSpackmanSet    SpackmanSetup( pgon, numverts, p_numrec )
  670. double    pgon[][2] ;
  671. int    numverts ;
  672. int    *p_numrec ;
  673. {
  674. int    p1, p2, degen ;
  675. double    denom, u1, v1, *pv[3] ;
  676. pSpackmanSet    pss, pss_return ;
  677. #ifdef    SORT
  678. double    u[2], v[2], len[2], *pv_temp ;
  679. #endif
  680.  
  681.     pss = pss_return =
  682.         (pSpackmanSet)malloc( (numverts-2) * sizeof( SpackmanSet )) ;
  683.     MALLOC_CHECK( pss ) ;
  684.  
  685.     degen = 0 ;
  686.  
  687.     for ( p1 = 1, p2 = 2 ; p2 < numverts ; p1++, p2++ ) {
  688.  
  689.     pv[0] = pgon[0] ;
  690.     pv[1] = pgon[p1] ;
  691.     pv[2] = pgon[p2] ;
  692.  
  693. #ifdef    SORT
  694.     /* Note that sorting can cause a mismatch of alpha/beta inequality
  695.      * tests.  In other words, test points on an interior line between
  696.      * test triangles will often then be wrong.
  697.      */
  698.     u[0] = pv[1][X] - pv[0][X] ;
  699.     u[1] = pv[2][X] - pv[0][X] ;
  700.     v[0] = pv[1][Y] - pv[0][Y] ;
  701.     v[1] = pv[2][Y] - pv[0][Y] ;
  702.     len[0] = u[0] * u[0] + v[0] * v[0] ;
  703.     len[1] = u[1] * u[1] + v[1] * v[1] ;
  704.  
  705.     /* compare two edges touching anchor point and put longest first */
  706.     /* we don't sort all three edges because the anchor point and
  707.      * values computed from it gets used for all triangles in the fan.
  708.      */
  709.     if ( len[0] < len[1] ) {
  710.         pv_temp = pv[1] ;
  711.         pv[1] = pv[2] ;
  712.         pv[2] = pv_temp ;
  713.     }
  714. #endif
  715.  
  716.     u1 = pv[1][X] - pv[0][X] ;
  717.     pss->u2 = pv[2][X] - pv[0][X] ;
  718.     v1 = pv[1][Y] - pv[0][Y] ;
  719.     pss->v2 = pv[2][Y] - pv[0][Y] ;
  720.     pss->u1_nonzero = !( u1 == 0.0 ) ;
  721.     if ( pss->u1_nonzero ) {
  722.         /* not zero, so compute inverse */
  723.         pss->inv_u1 = 1.0 / u1 ;
  724.         denom = pss->v2 * u1 - pss->u2 * v1 ;
  725.         if ( denom == 0.0 ) {
  726.         /* degenerate triangle, ignore it */
  727.         degen++ ;
  728.         goto Skip ;
  729.         } else {
  730.         pss->u1p = u1 / denom ;
  731.         pss->v1p = v1 / denom ;
  732.         }
  733.     } else {
  734.         if ( pss->u2 == 0.0 ) {
  735.         /* degenerate triangle, ignore it */
  736.         degen++ ;
  737.         goto Skip ;
  738.         } else {
  739.         /* not zero, so compute inverse */
  740.         pss->inv_u2 = 1.0 / pss->u2 ;
  741.         if ( v1 == 0.0 ) {
  742.             /* degenerate triangle, ignore it */
  743.             degen++ ;
  744.             goto Skip ;
  745.         } else {
  746.             pss->inv_v1 = 1.0 / v1 ;
  747.         }
  748.         }
  749.     }
  750.  
  751.     pss++ ;
  752.     Skip: ;
  753.     }
  754.  
  755.     /* number of Spackman records */
  756.     *p_numrec = numverts - degen - 2 ;
  757.     if ( degen ) {
  758.     pss = pss_return =
  759.         (pSpackmanSet)realloc( pss_return,
  760.             (numverts-2-degen) * sizeof( SpackmanSet )) ;
  761.     }
  762.  
  763.     return( pss_return ) ;
  764. }
  765.  
  766. /* barycentric, a la Gems I and Spackman's normalization precompute */
  767. int SpackmanTest( anchor, p_spackman_set, numrec, point )
  768. double    anchor[2] ;
  769. pSpackmanSet    p_spackman_set ;
  770. int    numrec ;
  771. double    point[2] ;
  772. {
  773. register pSpackmanSet    pss ;
  774. register int    inside_flag ;
  775. register int    nr ;
  776. register double tx, ty, vx0, vy0, u0, v0, alpha, beta ;
  777.  
  778.     tx = point[X] ;
  779.     ty = point[Y] ;
  780.     /* note that we really need only the first vertex of the polygon,
  781.      * so do not really need to keep the whole polygon around.
  782.      */
  783.     vx0 = anchor[X] ;
  784.     vy0 = anchor[Y] ;
  785.     u0 = tx - vx0 ;
  786.     v0 = ty - vy0 ;
  787.  
  788.     inside_flag = 0 ;
  789.  
  790.     for ( pss = p_spackman_set, nr = numrec+1 ; --nr ; pss++ ) {
  791.  
  792.     if ( pss->u1_nonzero ) {
  793.         /* 0 and 2 vertices have different X value */
  794.  
  795.         /* compute beta and check bounds */
  796.         /* we use "<= 0.0" so that points on the shared interior edge
  797.          * will (generally) be inside only one polygon.
  798.          */
  799.         beta = ( v0 * pss->u1p - u0 * pss->v1p ) ;
  800.         if ( ( beta <= 0.0 ) || ( beta > 1.0 ) ||
  801.  
  802.          /* compute alpha & check bounds */
  803.          ( ( alpha = ( u0 - beta * pss->u2 ) * pss->inv_u1 )
  804.             < 0.0 ) ) {
  805.  
  806.         /* whew! missed! */
  807.         goto NextTri ;
  808.         }
  809.     } else {
  810.         /* 0 and 2 vertices have same X value */
  811.  
  812.         /* compute beta and check bounds */
  813.         /* we use "<= 0.0" so that points on the shared interior edge
  814.          * will (generally) be inside only one polygon.
  815.          */
  816.         beta = u0 * pss->inv_u2 ;
  817.         if ( ( beta <= 0.0 ) || ( beta >= 1.0 ) ||
  818.  
  819.          /* compute alpha and check bounds */
  820.          ( ( alpha = ( v0 - beta * pss->v2 ) * pss->inv_v1 )
  821.             < 0.0 ) ) {
  822.  
  823.         /* whew! missed! */
  824.         goto NextTri ;
  825.         }
  826.     }
  827.  
  828.     /* check gamma */
  829.     if ( alpha + beta <= 1.0 ) {
  830.         /* survived */
  831.         inside_flag = !inside_flag ;
  832.     }
  833.  
  834.     NextTri: ;
  835.     }
  836.  
  837.     return( inside_flag ) ;
  838. }
  839.  
  840. void SpackmanCleanup( p_spackman_set )
  841. pSpackmanSet    p_spackman_set ;
  842. {
  843.     free( p_spackman_set ) ;
  844. }
  845.  
  846. /* ======= Trapezoid (bin) algorithm ======================================= */
  847.  
  848. /* Split polygons along set of y bins and sorts the edge fragments.  Testing
  849.  * is done against these fragments.
  850.  *
  851.  * Call setup with 2D polygon _pgon_ with _numverts_ number of vertices, the
  852.  * number of bins desired _bins_, and a pointer to a trapezoid structure
  853.  * _p_trap_set_.
  854.  * Call testing procedure with 2D polygon _pgon_ with _numverts_ number of
  855.  * vertices, _p_trap_set_ pointer to trapezoid structure, and test point
  856.  * _point_, returns 1 if inside, 0 if outside.
  857.  * Call cleanup with pointer to trapezoid structure to free space.
  858.  */
  859. void TrapezoidSetup( pgon, numverts, bins, p_trap_set )
  860. double    pgon[][2] ;
  861. int    numverts ;
  862. int    bins ;
  863. pTrapezoidSet    p_trap_set ;
  864. {
  865. double    *vtx0, *vtx1, *vtxa, *vtxb, slope ;
  866. int    i, j, bin_tot[TOT_BINS], ba, bb, id, full_cross, count ;
  867. double    fba, fbb, vx0, vx1, dy, vy0 ;
  868.  
  869.     p_trap_set->bins = bins ;
  870.     p_trap_set->trapz = (pTrapezoid)malloc( p_trap_set->bins *
  871.         sizeof(Trapezoid));
  872.     MALLOC_CHECK( p_trap_set->trapz ) ;
  873.  
  874.     p_trap_set->minx =
  875.     p_trap_set->maxx = pgon[0][X] ;
  876.     p_trap_set->miny =
  877.     p_trap_set->maxy = pgon[0][Y] ;
  878.  
  879.     for ( i = 1 ; i < numverts ; i++ ) {
  880.     if ( p_trap_set->minx > (vx0 = pgon[i][X]) ) {
  881.         p_trap_set->minx = vx0 ;
  882.     } else if ( p_trap_set->maxx < vx0 ) {
  883.         p_trap_set->maxx = vx0 ;
  884.     }
  885.  
  886.     if ( p_trap_set->miny > (vy0 = pgon[i][Y]) ) {
  887.         p_trap_set->miny = vy0 ;
  888.     } else if ( p_trap_set->maxy < vy0 ) {
  889.         p_trap_set->maxy = vy0 ;
  890.     }
  891.     }
  892.  
  893.     /* add a little to the bounds to ensure everything falls inside area */
  894.     p_trap_set->miny -= EPSILON * (p_trap_set->maxy-p_trap_set->miny) ;
  895.     p_trap_set->maxy += EPSILON * (p_trap_set->maxy-p_trap_set->miny) ;
  896.  
  897.     p_trap_set->ydelta =
  898.         (p_trap_set->maxy-p_trap_set->miny) / (double)p_trap_set->bins ;
  899.     p_trap_set->inv_ydelta = 1.0 / p_trap_set->ydelta ;
  900.  
  901.     /* find how many locations to allocate for each bin */
  902.     for ( i = 0 ; i < p_trap_set->bins ; i++ ) {
  903.     bin_tot[i] = 0 ;
  904.     }
  905.  
  906.     vtx0 = pgon[numverts-1] ;
  907.     for ( i = 0 ; i < numverts ; i++ ) {
  908.     vtx1 = pgon[i] ;
  909.  
  910.     /* skip if Y's identical (edge has no effect) */
  911.     if ( vtx0[Y] != vtx1[Y] ) {
  912.  
  913.         if ( vtx0[Y] < vtx1[Y] ) {
  914.         vtxa = vtx0 ;
  915.         vtxb = vtx1 ;
  916.         } else {
  917.         vtxa = vtx1 ;
  918.         vtxb = vtx0 ;
  919.         }
  920.         ba = (int)(( vtxa[Y]-p_trap_set->miny ) * p_trap_set->inv_ydelta) ;
  921.         fbb = ( vtxb[Y] - p_trap_set->miny ) * p_trap_set->inv_ydelta ;
  922.         bb = (int)fbb ;
  923.         /* if high vertex ends on a boundary, don't go into next boundary */
  924.         if ( fbb - (double)bb == 0.0 ) {
  925.         bb-- ;
  926.         }
  927.  
  928.         /* mark the bins with this edge */
  929.         for ( j = ba ; j <= bb ; j++ ) {
  930.         bin_tot[j]++ ;
  931.         }
  932.     }
  933.  
  934.     vtx0 = vtx1 ;
  935.     }
  936.  
  937.     /* allocate the bin contents and fill in some basics */
  938.     for ( i = 0 ; i < p_trap_set->bins ; i++ ) {
  939.     p_trap_set->trapz[i].edge_set =
  940.         (pEdge*)malloc( bin_tot[i] * sizeof(pEdge) ) ;
  941.     MALLOC_CHECK( p_trap_set->trapz[i].edge_set ) ;
  942.     for ( j = 0 ; j < bin_tot[i] ; j++ ) {
  943.         p_trap_set->trapz[i].edge_set[j] =
  944.         (pEdge)malloc( sizeof(Edge) ) ;
  945.         MALLOC_CHECK( p_trap_set->trapz[i].edge_set[j] ) ;
  946.     }
  947.  
  948.     /* start these off at some awful values; refined below */
  949.     p_trap_set->trapz[i].minx = p_trap_set->maxx ;
  950.     p_trap_set->trapz[i].maxx = p_trap_set->minx ;
  951.     p_trap_set->trapz[i].count = 0 ;
  952.     }
  953.  
  954.     /* now go through list yet again, putting edges in bins */
  955.     vtx0 = pgon[numverts-1] ;
  956.     id = numverts-1 ;
  957.     for ( i = 0 ; i < numverts ; i++ ) {
  958.     vtx1 = pgon[i] ;
  959.  
  960.     /* we can skip edge if Y's are equal */
  961.     if ( vtx0[Y] != vtx1[Y] ) {
  962.         if ( vtx0[Y] < vtx1[Y] ) {
  963.         vtxa = vtx0 ;
  964.         vtxb = vtx1 ;
  965.         } else {
  966.         vtxa = vtx1 ;
  967.         vtxb = vtx0 ;
  968.         }
  969.         fba = ( vtxa[Y] - p_trap_set->miny ) * p_trap_set->inv_ydelta ;
  970.         ba = (int)fba ;
  971.         fbb = ( vtxb[Y] - p_trap_set->miny ) * p_trap_set->inv_ydelta ;
  972.         bb = (int)fbb ;
  973.         /* if high vertex ends on a boundary, don't go into it */
  974.         if ( fbb == (double)bb ) {
  975.         bb-- ;
  976.         }
  977.  
  978.         vx0 = vtxa[X] ;
  979.         dy = vtxa[Y] - vtxb[Y] ;
  980.         slope = p_trap_set->ydelta * ( vtxa[X] - vtxb[X] ) / dy ;
  981.  
  982.         /* set vx1 in case loop is not entered */
  983.         vx1 = vx0 ;
  984.         full_cross = 0 ;
  985.  
  986.         for ( j = ba ; j < bb ; j++, vx0 = vx1 ) {
  987.         /* could increment vx1, but for greater accuracy recompute it */
  988.         vx1 = vtxa[X] + ( (double)(j+1) - fba ) * slope ;
  989.  
  990.         count = p_trap_set->trapz[j].count++ ;
  991.         p_trap_set->trapz[j].edge_set[count]->id = id ;
  992.         p_trap_set->trapz[j].edge_set[count]->full_cross = full_cross;
  993.         TrapBound( j, count, vx0, vx1, p_trap_set ) ;
  994.         full_cross = 1 ;
  995.         }
  996.  
  997.         /* at last bin - fill as above, but with vx1 = vtxb[X] */
  998.         vx0 = vx1 ;
  999.         vx1 = vtxb[X] ;
  1000.         count = p_trap_set->trapz[bb].count++ ;
  1001.         p_trap_set->trapz[bb].edge_set[count]->id = id ;
  1002.         /* the last bin is never a full crossing */
  1003.         p_trap_set->trapz[bb].edge_set[count]->full_cross = 0 ;
  1004.         TrapBound( bb, count, vx0, vx1, p_trap_set ) ;
  1005.     }
  1006.  
  1007.     vtx0 = vtx1 ;
  1008.     id = i ;
  1009.     }
  1010.  
  1011.     /* finally, sort the bins' contents by minx */
  1012.     for ( i = 0 ; i < p_trap_set->bins ; i++ ) {
  1013.     qsort( p_trap_set->trapz[i].edge_set, p_trap_set->trapz[i].count,
  1014.         sizeof(pEdge), CompareEdges ) ;
  1015.     }
  1016. }
  1017.  
  1018. void TrapBound( j, count, vx0, vx1, p_trap_set )
  1019. int    j, count ;
  1020. double    vx0, vx1 ;
  1021. pTrapezoidSet    p_trap_set ;
  1022. {
  1023. double    xt ;
  1024.  
  1025.     if ( vx0 > vx1 ) {
  1026.     xt = vx0 ;
  1027.     vx0 = vx1 ;
  1028.     vx1 = xt ;
  1029.     }
  1030.  
  1031.     if ( p_trap_set->trapz[j].minx > vx0 ) {
  1032.     p_trap_set->trapz[j].minx = vx0 ;
  1033.     }
  1034.     if ( p_trap_set->trapz[j].maxx < vx1 ) {
  1035.     p_trap_set->trapz[j].maxx = vx1 ;
  1036.     }
  1037.     p_trap_set->trapz[j].edge_set[count]->minx = vx0 ;
  1038.     p_trap_set->trapz[j].edge_set[count]->maxx = vx1 ;
  1039. }
  1040.  
  1041. /* used by qsort to sort */
  1042. int CompareEdges( u, v )
  1043. pEdge    *u, *v ;
  1044. {
  1045.     if ( (*u)->minx == (*v)->minx ) {
  1046.     return( 0 ) ;
  1047.     } else {
  1048.     return( (*u)->minx < (*v)->minx ? -1 : 1 ) ;
  1049.     }
  1050. }
  1051.  
  1052. int TrapezoidTest( pgon, numverts, p_trap_set, point )
  1053. double    pgon[][2] ;
  1054. int    numverts ;
  1055. pTrapezoidSet    p_trap_set ;
  1056. double    point[2] ;
  1057. {
  1058. int    j, b, count, id ;
  1059. double    tx, ty, *vtx0, *vtx1 ;
  1060. pEdge    *pp_bin ;
  1061. pTrapezoid    p_trap ;
  1062. int    inside_flag ;
  1063.  
  1064.     inside_flag = 0 ;
  1065.  
  1066.     /* first, is point inside bounding rectangle? */
  1067.     if ( ( ty = point[Y] ) < p_trap_set->miny ||
  1068.      ty >= p_trap_set->maxy ||
  1069.      ( tx = point[X] ) < p_trap_set->minx ||
  1070.      tx >= p_trap_set->maxx ) {
  1071.  
  1072.     /* outside of box */
  1073.     return( 0 ) ;
  1074.     }
  1075.  
  1076.     /* what bin are we in? */
  1077.     b = ( ty - p_trap_set->miny ) * p_trap_set->inv_ydelta ;
  1078.  
  1079.     /* find if we're inside this bin's bounds */
  1080.     if ( tx < (p_trap = &p_trap_set->trapz[b])->minx ||
  1081.      tx > p_trap->maxx ) {
  1082.  
  1083.     /* outside of box */
  1084.     return( 0 ) ;
  1085.     }
  1086.  
  1087.     /* now search bin for crossings */
  1088.     pp_bin = p_trap->edge_set ;
  1089.     count = p_trap->count ;
  1090.     for ( j = 0 ; j < count ; j++, pp_bin++ ) {
  1091.     if ( tx < (*pp_bin)->minx ) {
  1092.  
  1093.         /* all remaining edges are to right of point, so test them */
  1094.         do {
  1095.         if ( (*pp_bin)->full_cross ) {
  1096.             inside_flag = !inside_flag ;
  1097.         } else {
  1098.             id = (*pp_bin)->id ;
  1099.             if ( ( ty <= pgon[id][Y] ) !=
  1100.                 ( ty <= pgon[(id+1)%numverts][Y] ) ) {
  1101.  
  1102.             /* point crosses edge in Y, so must cross */
  1103.             inside_flag = !inside_flag ;
  1104.             }
  1105.         }
  1106.         pp_bin++ ;
  1107.         } while ( ++j < count ) ;
  1108.         goto Exit;
  1109.  
  1110.     } else if ( tx < (*pp_bin)->maxx ) {
  1111.         /* edge is overlapping point in X, check it */
  1112.         id = (*pp_bin)->id ;
  1113.         vtx0 = pgon[id] ;
  1114.         vtx1 = pgon[(id+1)%numverts] ;
  1115.  
  1116.         if ( (*pp_bin)->full_cross ||
  1117.          ( ty <= vtx0[Y] ) != ( ty <= vtx1[Y] ) ) {
  1118.  
  1119.         /* edge crosses in Y, so have to do full crossings test */
  1120.         if ( (vtx0[X] -
  1121.             (vtx0[Y] - ty )*
  1122.             ( vtx1[X]-vtx0[X])/(vtx1[Y]-vtx0[Y])) >= tx ) {
  1123.             inside_flag = !inside_flag ;
  1124.         }
  1125.         }
  1126.  
  1127.     } /* else edge is to left of point, ignore it */
  1128.     }
  1129.  
  1130.     Exit:
  1131.     return( inside_flag ) ;
  1132. }
  1133.  
  1134. void TrapezoidCleanup( p_trap_set )
  1135. pTrapezoidSet    p_trap_set ;
  1136. {
  1137. int    i, j, count ;
  1138.  
  1139.     for ( i = 0 ; i < p_trap_set->bins ; i++ ) {
  1140.     /* all of these should have bin sets, but check just in case */
  1141.     if ( p_trap_set->trapz[i].edge_set ) {
  1142.         count = p_trap_set->trapz[i].count ;
  1143.         for ( j = 0 ; j < count ; j++ ) {
  1144.         if ( p_trap_set->trapz[i].edge_set[j] ) {
  1145.             free( p_trap_set->trapz[i].edge_set[j] ) ;
  1146.         }
  1147.         }
  1148.         free( p_trap_set->trapz[i].edge_set ) ;
  1149.     }
  1150.     }
  1151.     free( p_trap_set->trapz ) ;
  1152. }
  1153.  
  1154. /* ======= Grid algorithm ================================================= */
  1155.  
  1156. /* Impose a grid upon the polygon and test only the local edges against the
  1157.  * point.
  1158.  *
  1159.  * Call setup with 2D polygon _pgon_ with _numverts_ number of vertices,
  1160.  * grid resolution _resolution_ and a pointer to a grid structure _p_gs_.
  1161.  * Call testing procedure with a pointer to this array and test point _point_,
  1162.  * returns 1 if inside, 0 if outside.
  1163.  * Call cleanup with pointer to grid structure to free space.
  1164.  */
  1165.  
  1166. /* Strategy for setup:
  1167.  *   Get bounds of polygon, allocate grid.
  1168.  *   "Walk" each edge of the polygon and note which edges have been crossed
  1169.  *     and what cells are entered (points on a grid edge are always considered
  1170.  *     to be above that edge).    Keep a record of the edges overlapping a cell.
  1171.  *     For cells with edges, determine if any cell border has no edges passing
  1172.  *     through it and so can be used for shooting a test ray.
  1173.  *     Keep track of the parity of the x (horizontal) grid cell borders for
  1174.  *     use in determining whether the grid corners are inside or outside.
  1175.  */
  1176. void GridSetup( pgon, numverts, resolution, p_gs )
  1177. double    pgon[][2] ;
  1178. int    numverts ;
  1179. int    resolution ;
  1180. pGridSet    p_gs ;
  1181. {
  1182. double    *vtx0, *vtx1, *vtxa, *vtxb, *p_gl ;
  1183. int    i, j, gc_clear_flags ;
  1184. double    vx0, vx1, vy0, vy1, gxdiff, gydiff, eps ;
  1185. pGridCell    p_gc, p_ngc ;
  1186. double    xdiff, ydiff, tmax, inv_x, inv_y, xdir, ydir, t_near, tx, ty ;
  1187. double    tgcx, tgcy ;
  1188. int    gcx, gcy, sign_x ;
  1189. int    y_flag, io_state ;
  1190.  
  1191.     p_gs->xres = p_gs->yres = resolution ;
  1192.     p_gs->tot_cells = p_gs->xres * p_gs->yres ;
  1193.     p_gs->glx = (double *)malloc( (p_gs->xres+1) * sizeof(double));
  1194.     MALLOC_CHECK( p_gs->glx ) ;
  1195.     p_gs->gly = (double *)malloc( (p_gs->yres+1) * sizeof(double));
  1196.     MALLOC_CHECK( p_gs->gly ) ;
  1197.     p_gs->gc = (pGridCell)malloc( p_gs->tot_cells * sizeof(GridCell));
  1198.     MALLOC_CHECK( p_gs->gc ) ;
  1199.  
  1200.     p_gs->minx =
  1201.     p_gs->maxx = pgon[0][X] ;
  1202.     p_gs->miny =
  1203.     p_gs->maxy = pgon[0][Y] ;
  1204.  
  1205.     /* find bounds of polygon */
  1206.     for ( i = 1 ; i < numverts ; i++ ) {
  1207.     vx0 = pgon[i][X] ;
  1208.     if ( p_gs->minx > vx0 ) {
  1209.         p_gs->minx = vx0 ;
  1210.     } else if ( p_gs->maxx < vx0 ) {
  1211.         p_gs->maxx = vx0 ;
  1212.     }
  1213.  
  1214.     vy0 = pgon[i][Y] ;
  1215.     if ( p_gs->miny > vy0 ) {
  1216.         p_gs->miny = vy0 ;
  1217.     } else if ( p_gs->maxy < vy0 ) {
  1218.         p_gs->maxy = vy0 ;
  1219.     }
  1220.     }
  1221.  
  1222.     /* add a little to the bounds to ensure everything falls inside area */
  1223.     gxdiff = p_gs->maxx - p_gs->minx ;
  1224.     gydiff = p_gs->maxy - p_gs->miny ;
  1225.     p_gs->minx -= EPSILON * gxdiff ;
  1226.     p_gs->maxx += EPSILON * gxdiff ;
  1227.     p_gs->miny -= EPSILON * gydiff ;
  1228.     p_gs->maxy += EPSILON * gydiff ;
  1229.  
  1230.     /* avoid roundoff problems near corners by not getting too close to them */
  1231.     eps = 1e-9 * ( gxdiff + gydiff ) ;
  1232.  
  1233.     /* use the new bounds to compute cell widths */
  1234.     TryAgain:
  1235.     p_gs->xdelta =
  1236.         (p_gs->maxx-p_gs->minx) / (double)p_gs->xres ;
  1237.     p_gs->inv_xdelta = 1.0 / p_gs->xdelta ;
  1238.  
  1239.     p_gs->ydelta =
  1240.         (p_gs->maxy-p_gs->miny) / (double)p_gs->yres ;
  1241.     p_gs->inv_ydelta = 1.0 / p_gs->ydelta ;
  1242.  
  1243.     for ( i = 0, p_gl = p_gs->glx ; i < p_gs->xres ; i++ ) {
  1244.     *p_gl++ = p_gs->minx + i * p_gs->xdelta ;
  1245.     }
  1246.     /* make last grid corner precisely correct */
  1247.     *p_gl = p_gs->maxx ;
  1248.  
  1249.     for ( i = 0, p_gl = p_gs->gly ; i < p_gs->yres ; i++ ) {
  1250.     *p_gl++ = p_gs->miny + i * p_gs->ydelta ;
  1251.     }
  1252.     *p_gl = p_gs->maxy ;
  1253.  
  1254.     for ( i = 0, p_gc = p_gs->gc ; i < p_gs->tot_cells ; i++, p_gc++ ) {
  1255.     p_gc->tot_edges = 0 ;
  1256.     p_gc->gc_flags = 0x0 ;
  1257.     p_gc->gr = NULL ;
  1258.     }
  1259.  
  1260.     /* loop through edges and insert into grid structure */
  1261.     vtx0 = pgon[numverts-1] ;
  1262.     for ( i = 0 ; i < numverts ; i++ ) {
  1263.     vtx1 = pgon[i] ;
  1264.  
  1265.     if ( vtx0[Y] < vtx1[Y] ) {
  1266.         vtxa = vtx0 ;
  1267.         vtxb = vtx1 ;
  1268.     } else {
  1269.         vtxa = vtx1 ;
  1270.         vtxb = vtx0 ;
  1271.     }
  1272.  
  1273.     /* Set x variable for the direction of the ray */
  1274.     xdiff = vtxb[X] - vtxa[X] ;
  1275.     ydiff = vtxb[Y] - vtxa[Y] ;
  1276.     tmax = sqrt( xdiff * xdiff + ydiff * ydiff ) ;
  1277.  
  1278.     /* if edge is of 0 length, ignore it (useless edge) */
  1279.     if ( tmax == 0.0 ) goto NextEdge ;
  1280.  
  1281.     inv_y = tmax / ydiff ;
  1282.     xdir = xdiff / tmax ;
  1283.     ydir = ydiff / tmax ;
  1284.  
  1285.     gcx = (int)(( vtxa[X] - p_gs->minx ) * p_gs->inv_xdelta) ;
  1286.     gcy = (int)(( vtxa[Y] - p_gs->miny ) * p_gs->inv_ydelta) ;
  1287.  
  1288.     /* get information about slopes of edge, etc */
  1289.     if ( vtxa[X] == vtxb[X] ) {
  1290.         sign_x = 0 ;
  1291.         tx = HUGE ;
  1292.     } else {
  1293.         inv_x = tmax / xdiff ;
  1294.         tx = p_gs->xdelta * (double)gcx + p_gs->minx - vtxa[X] ;
  1295.         if ( vtxa[X] < vtxb[X] ) {
  1296.         sign_x = 1 ;
  1297.         tx += p_gs->xdelta ;
  1298.         tgcx = p_gs->xdelta * inv_x ;
  1299.         } else {
  1300.         sign_x = -1 ;
  1301.         tgcx = -p_gs->xdelta * inv_x ;
  1302.         }
  1303.         tx *= inv_x ;
  1304.     }
  1305.  
  1306.     if ( vtxa[Y] == vtxb[Y] ) {
  1307.         ty = HUGE ;
  1308.     } else {
  1309.         inv_y = tmax / ydiff ;
  1310.         ty = (p_gs->ydelta * (double)(gcy+1) + p_gs->miny - vtxa[Y])
  1311.         * inv_y ;
  1312.         tgcy = p_gs->ydelta * inv_y ;
  1313.     }
  1314.  
  1315.     p_gc = &p_gs->gc[gcy*p_gs->xres+gcx] ;
  1316.  
  1317.     vx0 = vtxa[X] ;
  1318.     vy0 = vtxa[Y] ;
  1319.  
  1320.     t_near = 0.0 ;
  1321.  
  1322.     do {
  1323.         /* choose the next boundary, but don't move yet */
  1324.         if ( tx <= ty ) {
  1325.         gcx += sign_x ;
  1326.  
  1327.         ty -= tx ;
  1328.         t_near += tx ;
  1329.         tx = tgcx ;
  1330.  
  1331.         /* note which edge is hit when leaving this cell */
  1332.         if ( t_near < tmax ) {
  1333.             if ( sign_x > 0 ) {
  1334.             p_gc->gc_flags |= GC_R_EDGE_HIT ;
  1335.             vx1 = p_gs->glx[gcx] ;
  1336.             } else {
  1337.             p_gc->gc_flags |= GC_L_EDGE_HIT ;
  1338.             vx1 = p_gs->glx[gcx+1] ;
  1339.             }
  1340.  
  1341.             /* get new location */
  1342.             vy1 = t_near * ydir + vtxa[Y] ;
  1343.         } else {
  1344.             /* end of edge, so get exact value */
  1345.             vx1 = vtxb[X] ;
  1346.             vy1 = vtxb[Y] ;
  1347.         }
  1348.  
  1349.         y_flag = FALSE ;
  1350.  
  1351.         } else {
  1352.  
  1353.         gcy++ ;
  1354.  
  1355.         tx -= ty ;
  1356.         t_near += ty ;
  1357.         ty = tgcy ;
  1358.  
  1359.         /* note top edge is hit when leaving this cell */
  1360.         if ( t_near < tmax ) {
  1361.             p_gc->gc_flags |= GC_T_EDGE_HIT ;
  1362.             /* this toggles the parity bit */
  1363.             p_gc->gc_flags ^= GC_T_EDGE_PARITY ;
  1364.  
  1365.             /* get new location */
  1366.             vx1 = t_near * xdir + vtxa[X] ;
  1367.             vy1 = p_gs->gly[gcy] ;
  1368.         } else {
  1369.             /* end of edge, so get exact value */
  1370.             vx1 = vtxb[X] ;
  1371.             vy1 = vtxb[Y] ;
  1372.         }
  1373.  
  1374.         y_flag = TRUE ;
  1375.         }
  1376.  
  1377.         /* check for corner crossing, then mark the cell we're in */
  1378.         if ( !AddGridRecAlloc( p_gc, vx0, vy0, vx1, vy1, eps ) ) {
  1379.         /* warning, danger - we have just crossed a corner.
  1380.          * There are all kinds of topological messiness we could
  1381.          * do to get around this case, but they're a headache.
  1382.          * The simplest recovery is just to change the extents a bit
  1383.          * and redo the meshing, so that hopefully no edges will
  1384.          * perfectly cross a corner.  Since it's a preprocess, we
  1385.          * don't care too much about the time to do it.
  1386.          */
  1387.  
  1388.         /* clean out all grid records */
  1389.         for ( i = 0, p_gc = p_gs->gc
  1390.             ; i < p_gs->tot_cells
  1391.             ; i++, p_gc++ ) {
  1392.  
  1393.             if ( p_gc->gr ) {
  1394.             free( p_gc->gr ) ;
  1395.             }
  1396.         }
  1397.  
  1398.         /* make the bounding box ever so slightly larger, hopefully
  1399.          * changing the alignment of the corners.
  1400.          */
  1401.         p_gs->minx -= EPSILON * gxdiff * 0.24 ;
  1402.         p_gs->miny -= EPSILON * gydiff * 0.10 ;
  1403.  
  1404.         /* yes, it's the dreaded goto - run in fear for your lives! */
  1405.         goto TryAgain ;
  1406.         }
  1407.  
  1408.         if ( t_near < tmax ) {
  1409.         /* note how we're entering the next cell */
  1410.         /* TBD: could be done faster by incrementing index in the
  1411.          * incrementing code, above */
  1412.         p_gc = &p_gs->gc[gcy*p_gs->xres+gcx] ;
  1413.  
  1414.         if ( y_flag ) {
  1415.             p_gc->gc_flags |= GC_B_EDGE_HIT ;
  1416.             /* this toggles the parity bit */
  1417.             p_gc->gc_flags ^= GC_B_EDGE_PARITY ;
  1418.         } else {
  1419.             p_gc->gc_flags |=
  1420.             ( sign_x > 0 ) ? GC_L_EDGE_HIT : GC_R_EDGE_HIT ;
  1421.         }
  1422.         }
  1423.  
  1424.         vx0 = vx1 ;
  1425.         vy0 = vy1 ;
  1426.     }
  1427.     /* have we gone further than the end of the edge? */
  1428.     while ( t_near < tmax ) ;
  1429.  
  1430.     NextEdge:
  1431.     vtx0 = vtx1 ;
  1432.     }
  1433.  
  1434.     /* the grid is all set up, now set up the inside/outside value of each
  1435.      * corner.
  1436.      */
  1437.     p_gc = p_gs->gc ;
  1438.     p_ngc = &p_gs->gc[p_gs->xres] ;
  1439.  
  1440.     /* we know the bottom and top rows are all outside, so no flag is set */
  1441.     for ( i = 1; i < p_gs->yres ; i++ ) {
  1442.     /* start outside */
  1443.     io_state = 0x0 ;
  1444.  
  1445.     for ( j = 0; j < p_gs->xres ; j++ ) {
  1446.  
  1447.         if ( io_state ) {
  1448.         /* change cell left corners to inside */
  1449.         p_gc->gc_flags |= GC_TL_IN ;
  1450.         p_ngc->gc_flags |= GC_BL_IN ;
  1451.         }
  1452.  
  1453.         if ( p_gc->gc_flags & GC_T_EDGE_PARITY ) {
  1454.         io_state = !io_state ;
  1455.         }
  1456.  
  1457.         if ( io_state ) {
  1458.         /* change cell right corners to inside */
  1459.         p_gc->gc_flags |= GC_TR_IN ;
  1460.         p_ngc->gc_flags |= GC_BR_IN ;
  1461.         }
  1462.  
  1463.         p_gc++ ;
  1464.         p_ngc++ ;
  1465.     }
  1466.     }
  1467.  
  1468.     p_gc = p_gs->gc ;
  1469.     for ( i = 0; i < p_gs->tot_cells ; i++ ) {
  1470.  
  1471.     /* reverse parity of edge clear (1==edge clear) */
  1472.     gc_clear_flags = p_gc->gc_flags ^ GC_ALL_EDGE_CLEAR ;
  1473.     if ( gc_clear_flags & GC_L_EDGE_CLEAR ) {
  1474.         p_gc->gc_flags |= GC_AIM_L ;
  1475.     } else
  1476.     if ( gc_clear_flags & GC_B_EDGE_CLEAR ) {
  1477.         p_gc->gc_flags |= GC_AIM_B ;
  1478.     } else
  1479.     if ( gc_clear_flags & GC_R_EDGE_CLEAR ) {
  1480.         p_gc->gc_flags |= GC_AIM_R ;
  1481.     } else
  1482.     if ( gc_clear_flags & GC_T_EDGE_CLEAR ) {
  1483.         p_gc->gc_flags |= GC_AIM_T ;
  1484.     } else {
  1485.         /* all edges have something on them, do full test */
  1486.         p_gc->gc_flags |= GC_AIM_C ;
  1487.     }
  1488.     p_gc++ ;
  1489.     }
  1490. }
  1491.  
  1492. int AddGridRecAlloc( p_gc, xa, ya, xb, yb, eps )
  1493. pGridCell    p_gc ;
  1494. double    xa,ya,xb,yb,eps ;
  1495. {
  1496. pGridRec    p_gr ;
  1497. double        slope, inv_slope ;
  1498.  
  1499.     if ( Near(ya, yb, eps) ) {
  1500.     if ( Near(xa, xb, eps) ) {
  1501.         /* edge is 0 length, so get rid of it */
  1502.         return( FALSE ) ;
  1503.     } else {
  1504.         /* horizontal line */
  1505.         slope = HUGE ;
  1506.         inv_slope = 0.0 ;
  1507.     }
  1508.     } else {
  1509.     if ( Near(xa, xb, eps) ) {
  1510.         /* vertical line */
  1511.         slope = 0.0 ;
  1512.         inv_slope = HUGE ;
  1513.     } else {
  1514.         slope = (xb-xa)/(yb-ya) ;
  1515.         inv_slope = (yb-ya)/(xb-xa) ;
  1516.     }
  1517.     }
  1518.  
  1519.     p_gc->tot_edges++ ;
  1520.     if ( p_gc->tot_edges <= 1 ) {
  1521.     p_gc->gr = (pGridRec)malloc( sizeof(GridRec) ) ;
  1522.     } else {
  1523.     p_gc->gr = (pGridRec)realloc( p_gc->gr,
  1524.         p_gc->tot_edges * sizeof(GridRec) ) ;
  1525.     }
  1526.     MALLOC_CHECK( p_gc->gr ) ;
  1527.     p_gr = &p_gc->gr[p_gc->tot_edges-1] ;
  1528.  
  1529.     p_gr->slope = slope ;
  1530.     p_gr->inv_slope = inv_slope ;
  1531.  
  1532.     p_gr->xa = xa ;
  1533.     p_gr->ya = ya ;
  1534.     if ( xa <= xb ) {
  1535.     p_gr->minx = xa ;
  1536.     p_gr->maxx = xb ;
  1537.     } else {
  1538.     p_gr->minx = xb ;
  1539.     p_gr->maxx = xa ;
  1540.     }
  1541.     if ( ya <= yb ) {
  1542.     p_gr->miny = ya ;
  1543.     p_gr->maxy = yb ;
  1544.     } else {
  1545.     p_gr->miny = yb ;
  1546.     p_gr->maxy = ya ;
  1547.     }
  1548.  
  1549.     /* P2 - P1 */
  1550.     p_gr->ax = xb - xa ;
  1551.     p_gr->ay = yb - ya ;
  1552.  
  1553.     return( TRUE ) ;
  1554. }
  1555.  
  1556. /* Test point against grid and edges in the cell (if any).  Algorithm:
  1557.  *    Check bounding box; if outside then return.
  1558.  *    Check cell point is inside; if simple inside or outside then return.
  1559.  *    Find which edge or corner is considered to be the best for testing and
  1560.  *      send a test ray towards it, counting the crossings.  Add in the
  1561.  *      state of the edge or corner the ray went to and so determine the
  1562.  *      state of the point (inside or outside).
  1563.  */
  1564. int GridTest( p_gs, point )
  1565. register pGridSet    p_gs ;
  1566. double    point[2] ;
  1567. {
  1568. int    j, count, init_flag ;
  1569. pGridCell    p_gc ;
  1570. pGridRec    p_gr ;
  1571. double    tx, ty, xcell, ycell, bx,by,cx,cy, cornerx, cornery ;
  1572. double    alpha, beta, denom ;
  1573. unsigned short    gc_flags ;
  1574. int    inside_flag ;
  1575.  
  1576.     /* first, is point inside bounding rectangle? */
  1577.     if ( ( ty = point[Y] ) < p_gs->miny ||
  1578.      ty >= p_gs->maxy ||
  1579.      ( tx = point[X] ) < p_gs->minx ||
  1580.      tx >= p_gs->maxx ) {
  1581.  
  1582.     /* outside of box */
  1583.     inside_flag = FALSE ;
  1584.     } else {
  1585.  
  1586.     /* what cell are we in? */
  1587.     ycell = ( ty - p_gs->miny ) * p_gs->inv_ydelta ;
  1588.     xcell = ( tx - p_gs->minx ) * p_gs->inv_xdelta ;
  1589.     p_gc = &p_gs->gc[((int)ycell)*p_gs->xres + (int)xcell] ;
  1590.  
  1591.     /* is cell simple? */
  1592.     count = p_gc->tot_edges ;
  1593.     if ( count ) {
  1594.         /* no, so find an edge which is free. */
  1595.         gc_flags = p_gc->gc_flags ;
  1596.         p_gr = p_gc->gr ;
  1597.         switch( gc_flags & GC_AIM ) {
  1598.         case GC_AIM_L:
  1599.         /* left edge is clear, shoot X- ray */
  1600.         /* note - this next statement requires that GC_BL_IN is 1 */
  1601.         inside_flag = gc_flags & GC_BL_IN ;
  1602.         for ( j = count+1 ; --j ; p_gr++ ) {
  1603.             /* test if y is between edges */
  1604.             if ( ty >= p_gr->miny && ty < p_gr->maxy ) {
  1605.             if ( tx > p_gr->maxx ) {
  1606.                 inside_flag = !inside_flag ;
  1607.             } else if ( tx > p_gr->minx ) {
  1608.                 /* full computation */
  1609.                 if ( ( p_gr->xa -
  1610.                 ( p_gr->ya - ty ) * p_gr->slope ) < tx ) {
  1611.                 inside_flag = !inside_flag ;
  1612.                 }
  1613.             }
  1614.             }
  1615.         }
  1616.         break ;
  1617.  
  1618.         case GC_AIM_B:
  1619.         /* bottom edge is clear, shoot Y+ ray */
  1620.         /* note - this next statement requires that GC_BL_IN is 1 */
  1621.         inside_flag = gc_flags & GC_BL_IN ;
  1622.         for ( j = count+1 ; --j ; p_gr++ ) {
  1623.             /* test if x is between edges */
  1624.             if ( tx >= p_gr->minx && tx < p_gr->maxx ) {
  1625.             if ( ty > p_gr->maxy ) {
  1626.                 inside_flag = !inside_flag ;
  1627.             } else if ( ty > p_gr->miny ) {
  1628.                 /* full computation */
  1629.                 if ( ( p_gr->ya - ( p_gr->xa - tx ) *
  1630.                     p_gr->inv_slope ) < ty ) {
  1631.                 inside_flag = !inside_flag ;
  1632.                 }
  1633.             }
  1634.             }
  1635.         }
  1636.         break ;
  1637.  
  1638.         case GC_AIM_R:
  1639.         /* right edge is clear, shoot X+ ray */
  1640.         inside_flag = (gc_flags & GC_TR_IN) ? 1 : 0 ;
  1641.  
  1642.         /* TBD: Note, we could have sorted the edges to be tested
  1643.          * by miny or somesuch, and so be able to cut testing
  1644.          * short when the list's miny > point.y .
  1645.          */
  1646.         for ( j = count+1 ; --j ; p_gr++ ) {
  1647.             /* test if y is between edges */
  1648.             if ( ty >= p_gr->miny && ty < p_gr->maxy ) {
  1649.             if ( tx <= p_gr->minx ) {
  1650.                 inside_flag = !inside_flag ;
  1651.             } else if ( tx <= p_gr->maxx ) {
  1652.                 /* full computation */
  1653.                 if ( ( p_gr->xa -
  1654.                 ( p_gr->ya - ty ) * p_gr->slope ) >= tx ) {
  1655.                 inside_flag = !inside_flag ;
  1656.                 }
  1657.             }
  1658.             }
  1659.         }
  1660.         break ;
  1661.  
  1662.         case GC_AIM_T:
  1663.         /* top edge is clear, shoot Y+ ray */
  1664.         inside_flag = (gc_flags & GC_TR_IN) ? 1 : 0 ;
  1665.         for ( j = count+1 ; --j ; p_gr++ ) {
  1666.             /* test if x is between edges */
  1667.             if ( tx >= p_gr->minx && tx < p_gr->maxx ) {
  1668.             if ( ty <= p_gr->miny ) {
  1669.                 inside_flag = !inside_flag ;
  1670.             } else if ( ty <= p_gr->maxy ) {
  1671.                 /* full computation */
  1672.                 if ( ( p_gr->ya - ( p_gr->xa - tx ) *
  1673.                     p_gr->inv_slope ) >= ty ) {
  1674.                 inside_flag = !inside_flag ;
  1675.                 }
  1676.             }
  1677.             }
  1678.         }
  1679.         break ;
  1680.  
  1681.         case GC_AIM_C:
  1682.         /* no edge is clear, bite the bullet and test
  1683.          * against the bottom left corner.
  1684.          * We use Franklin Antonio's algorithm (Graphics Gems III).
  1685.          */
  1686.         /* TBD: Faster yet might be to test against the closest
  1687.          * corner to the cell location, but our hope is that we
  1688.          * rarely need to do this testing at all.
  1689.          */
  1690.         inside_flag = ((gc_flags & GC_BL_IN) == GC_BL_IN) ;
  1691.         init_flag = TRUE ;
  1692.  
  1693.         /* get lower left corner coordinate */
  1694.         cornerx = p_gs->glx[(int)xcell] ;
  1695.         cornery = p_gs->gly[(int)ycell] ;
  1696.         for ( j = count+1 ; --j ; p_gr++ ) {
  1697.  
  1698.             /* quick out test: if test point is
  1699.              * less than minx & miny, edge cannot overlap.
  1700.              */
  1701.             if ( tx >= p_gr->minx && ty >= p_gr->miny ) {
  1702.  
  1703.             /* quick test failed, now check if test point and
  1704.              * corner are on different sides of edge.
  1705.              */
  1706.             if ( init_flag ) {
  1707.                 /* Compute these at most once for test */
  1708.                 /* P3 - P4 */
  1709.                 bx = tx - cornerx ;
  1710.                 by = ty - cornery ;
  1711.                 init_flag = FALSE ;
  1712.             }
  1713.             denom = p_gr->ay * bx - p_gr->ax * by ;
  1714.             if ( denom != 0.0 ) {
  1715.                 /* lines are not collinear, so continue */
  1716.                 /* P1 - P3 */
  1717.                 cx = p_gr->xa - tx ;
  1718.                 cy = p_gr->ya - ty ;
  1719.                 alpha = by * cx - bx * cy ;
  1720.                 if ( denom > 0.0 ) {
  1721.                 if ( alpha < 0.0 || alpha >= denom ) {
  1722.                     /* test edge not hit */
  1723.                     goto NextEdge ;
  1724.                 }
  1725.                 beta = p_gr->ax * cy - p_gr->ay * cx ;
  1726.                 if ( beta < 0.0 || beta >= denom ) {
  1727.                     /* polygon edge not hit */
  1728.                     goto NextEdge ;
  1729.                 }
  1730.                 } else {
  1731.                 if ( alpha > 0.0 || alpha <= denom ) {
  1732.                     /* test edge not hit */
  1733.                     goto NextEdge ;
  1734.                 }
  1735.                 beta = p_gr->ax * cy - p_gr->ay * cx ;
  1736.                 if ( beta > 0.0 || beta <= denom ) {
  1737.                     /* polygon edge not hit */
  1738.                     goto NextEdge ;
  1739.                 }
  1740.                 }
  1741.                 inside_flag = !inside_flag ;
  1742.             }
  1743.  
  1744.             }
  1745.             NextEdge: ;
  1746.         }
  1747.         break ;
  1748.         }
  1749.     } else {
  1750.         /* simple cell, so if lower left corner is in,
  1751.          * then cell is inside.
  1752.          */
  1753.         inside_flag = p_gc->gc_flags & GC_BL_IN ;
  1754.     }
  1755.     }
  1756.  
  1757.     return( inside_flag ) ;
  1758. }
  1759.  
  1760. void GridCleanup( p_gs )
  1761. pGridSet    p_gs ;
  1762. {
  1763. int    i ;
  1764. pGridCell    p_gc ;
  1765.  
  1766.     for ( i = 0, p_gc = p_gs->gc
  1767.     ; i < p_gs->tot_cells
  1768.     ; i++, p_gc++ ) {
  1769.  
  1770.     if ( p_gc->gr ) {
  1771.         free( p_gc->gr ) ;
  1772.     }
  1773.     }
  1774.     free( p_gs->glx ) ;
  1775.     free( p_gs->gly ) ;
  1776.     free( p_gs->gc ) ;
  1777. }
  1778.  
  1779. /* ======= Exterior (convex only) algorithm =============================== */
  1780.  
  1781. /* Test the edges of the convex polygon against the point.  If the point is
  1782.  * outside any edge, the point is outside the polygon.
  1783.  *
  1784.  * Call setup with 2D polygon _pgon_ with _numverts_ number of vertices,
  1785.  * which returns a pointer to a plane set array.
  1786.  * Call testing procedure with a pointer to this array, _numverts_, and
  1787.  * test point _point_, returns 1 if inside, 0 if outside.
  1788.  * Call cleanup with pointer to plane set array to free space.
  1789.  *
  1790.  * RANDOM can be defined for this test.
  1791.  * CONVEX must be defined for this test; it is not usable for general polygons.
  1792.  */
  1793.  
  1794. #ifdef    CONVEX
  1795. /* make exterior plane set */
  1796. pPlaneSet ExteriorSetup( pgon, numverts )
  1797. double    pgon[][2] ;
  1798. int    numverts ;
  1799. {
  1800. int    p1, p2, flip_edge ;
  1801. pPlaneSet    pps, pps_return ;
  1802. #ifdef    RANDOM
  1803. int    i, ind ;
  1804. PlaneSet    ps_temp ;
  1805. #endif
  1806.  
  1807.     pps = pps_return =
  1808.         (pPlaneSet)malloc( numverts * sizeof( PlaneSet )) ;
  1809.     MALLOC_CHECK( pps ) ;
  1810.  
  1811.     /* take cross product of vertex to find handedness */
  1812.     flip_edge = (pgon[0][X] - pgon[1][X]) * (pgon[1][Y] - pgon[2][Y] ) >
  1813.         (pgon[0][Y] - pgon[1][Y]) * (pgon[1][X] - pgon[2][X] ) ;
  1814.  
  1815.     /* Generate half-plane boundary equations now for faster testing later.
  1816.      * vx & vy are the edge's normal, c is the offset from the origin.
  1817.      */
  1818.     for ( p1 = numverts-1, p2 = 0 ; p2 < numverts ; p1 = p2, p2++, pps++ ) {
  1819.     pps->vx = pgon[p1][Y] - pgon[p2][Y] ;
  1820.     pps->vy = pgon[p2][X] - pgon[p1][X] ;
  1821.     pps->c = pps->vx * pgon[p1][X] + pps->vy * pgon[p1][Y] ;
  1822.  
  1823.     /* check sense and reverse plane edge if need be */
  1824.     if ( flip_edge ) {
  1825.         pps->vx = -pps->vx ;
  1826.         pps->vy = -pps->vy ;
  1827.         pps->c  = -pps->c ;
  1828.     }
  1829.     }
  1830.  
  1831. #ifdef    RANDOM
  1832.     /* Randomize the order of the edges to improve chance of early out */
  1833.     /* There are better orders, but the default order is the worst */
  1834.     for ( i = 0, pps = pps_return
  1835.     ; i < numverts
  1836.     ; i++ ) {
  1837.  
  1838.     ind = (int)(RAN01() * numverts ) ;
  1839.     if ( ( ind < 0 ) || ( ind >= numverts ) ) {
  1840.         fprintf( stderr,
  1841.         "Yikes, the random number generator is returning values\n" ) ;
  1842.         fprintf( stderr,
  1843.         "outside the range [0.0,1.0), so please fix the code!\n" ) ;
  1844.         ind = 0 ;
  1845.     }
  1846.  
  1847.     /* swap edges */
  1848.     ps_temp = *pps ;
  1849.     *pps = pps_return[ind] ;
  1850.     pps_return[ind] = ps_temp ;
  1851.     }
  1852. #endif
  1853.  
  1854.     return( pps_return ) ;
  1855. }
  1856.  
  1857. /* Check point for outside of all planes */
  1858. /* note that we don't need "pgon", since it's been processed into
  1859.  * its corresponding PlaneSet.
  1860.  */
  1861. int ExteriorTest( p_ext_set, numverts, point )
  1862. pPlaneSet    p_ext_set ;
  1863. int    numverts ;
  1864. double    point[2] ;
  1865. {
  1866. register PlaneSet    *pps ;
  1867. register int    p0 ;
  1868. register double tx, ty ;
  1869. int    inside_flag ;
  1870.  
  1871.     tx = point[X] ;
  1872.     ty = point[Y] ;
  1873.  
  1874.     for ( p0 = numverts+1, pps = p_ext_set ; --p0 ; pps++ ) {
  1875.  
  1876.     /* test if the point is outside this edge */
  1877.     if ( pps->vx * tx + pps->vy * ty > pps->c ) {
  1878.         return( 0 ) ;
  1879.     }
  1880.     }
  1881.     /* if we make it to here, we were inside all edges */
  1882.     return( 1 ) ;
  1883. }
  1884.  
  1885. void ExteriorCleanup( p_ext_set )
  1886. pPlaneSet    p_ext_set ;
  1887. {
  1888.     free( p_ext_set ) ;
  1889. }
  1890. #endif
  1891.  
  1892. /* ======= Inclusion (convex only) algorithm ============================== */
  1893.  
  1894. /* Create an efficiency structure (see Preparata) for the convex polygon which
  1895.  * allows binary searching to find which edge to test the point against.  This
  1896.  * algorithm is O(log n).
  1897.  *
  1898.  * Call setup with 2D polygon _pgon_ with _numverts_ number of vertices,
  1899.  * which returns a pointer to an inclusion anchor structure.
  1900.  * Call testing procedure with a pointer to this structure and test point
  1901.  * _point_, returns 1 if inside, 0 if outside.
  1902.  * Call cleanup with pointer to inclusion anchor structure to free space.
  1903.  *
  1904.  * CONVEX must be defined for this test; it is not usable for general polygons.
  1905.  */
  1906.  
  1907. #ifdef    CONVEX
  1908. /* make inclusion wedge set */
  1909. pInclusionAnchor InclusionSetup( pgon, numverts )
  1910. double    pgon[][2] ;
  1911. int    numverts ;
  1912. {
  1913. int    pc, p1, p2, flip_edge ;
  1914. double    ax,ay, qx,qy, wx,wy, len ;
  1915. pInclusionAnchor pia ;
  1916. pInclusionSet    pis ;
  1917.  
  1918.     /* double the first edge to avoid needing modulo during test search */
  1919.     pia = (pInclusionAnchor)malloc( sizeof( InclusionAnchor )) ;
  1920.     MALLOC_CHECK( pia ) ;
  1921.     pis = pia->pis =
  1922.         (pInclusionSet)malloc( (numverts+1) * sizeof( InclusionSet )) ;
  1923.     MALLOC_CHECK( pis ) ;
  1924.  
  1925.     pia->hi_start = numverts - 1 ;
  1926.  
  1927.     /* get average point to make wedges from */
  1928.     qx = qy = 0.0 ;
  1929.     for ( p2 = 0 ; p2 < numverts ; p2++ ) {
  1930.     qx += pgon[p2][X] ;
  1931.     qy += pgon[p2][Y] ;
  1932.     }
  1933.     pia->qx = qx /= (double)numverts ;
  1934.     pia->qy = qy /= (double)numverts ;
  1935.  
  1936.     /* take cross product of vertex to find handedness */
  1937.     pia->flip_edge = flip_edge =
  1938.         (pgon[0][X] - pgon[1][X]) * (pgon[1][Y] - pgon[2][Y] ) >
  1939.         (pgon[0][Y] - pgon[1][Y]) * (pgon[1][X] - pgon[2][X] ) ;
  1940.  
  1941.  
  1942.     ax = pgon[0][X] - qx ;
  1943.     ay = pgon[0][Y] - qy ;
  1944.     len = sqrt( ax * ax + ay * ay ) ;
  1945.     if ( len == 0.0 ) {
  1946.     fprintf( stderr, "sorry, polygon for inclusion test is defective\n" ) ;
  1947.     exit(1) ;
  1948.     }
  1949.     pia->ax = ax /= len ;
  1950.     pia->ay = ay /= len ;
  1951.  
  1952.     /* loop through edges, and double last edge */
  1953.     for ( pc = p1 = 0, p2 = 1
  1954.     ; pc <= numverts
  1955.     ; pc++, p1 = p2, p2 = (++p2)%numverts, pis++ ) {
  1956.  
  1957.     /* wedge border */
  1958.     wx = pgon[p1][X] - qx ;
  1959.     wy = pgon[p1][Y] - qy ;
  1960.     len = sqrt( wx * wx + wy * wy ) ;
  1961.     wx /= len ;
  1962.     wy /= len ;
  1963.  
  1964.     /* cosine of angle from anchor border to wedge border */
  1965.     pis->dot = ax * wx + ay * wy ;
  1966.     /* sign from cross product */
  1967.     if ( ( ax * wy > ay * wx ) == flip_edge ) {
  1968.         pis->dot = -2.0 - pis->dot ;
  1969.     }
  1970.  
  1971.     /* edge */
  1972.     pis->ex = pgon[p1][Y] - pgon[p2][Y] ;
  1973.     pis->ey = pgon[p2][X] - pgon[p1][X] ;
  1974.     pis->ec = pis->ex * pgon[p1][X] + pis->ey * pgon[p1][Y] ;
  1975.  
  1976.     /* check sense and reverse plane eqns if need be */
  1977.     if ( flip_edge ) {
  1978.         pis->ex = -pis->ex ;
  1979.         pis->ey = -pis->ey ;
  1980.         pis->ec = -pis->ec ;
  1981.     }
  1982.     }
  1983.     /* set first angle a little > 1.0 and last < -3.0 just to be safe. */
  1984.     pia->pis[0].dot = -3.001 ;
  1985.     pia->pis[numverts].dot = 1.001 ;
  1986.  
  1987.     return( pia ) ;
  1988. }
  1989.  
  1990. /* Find wedge point is in by binary search, then test wedge */
  1991. int InclusionTest( pia, point )
  1992. pInclusionAnchor    pia ;
  1993. double    point[2] ;
  1994. {
  1995. register double tx, ty, len, dot ;
  1996. int    inside_flag, lo, hi, ind ;
  1997. pInclusionSet    pis ;
  1998.  
  1999.     tx = point[X] - pia->qx ;
  2000.     ty = point[Y] - pia->qy ;
  2001.     len = sqrt( tx * tx + ty * ty ) ;
  2002.     /* check if point is exactly at anchor point (which is inside polygon) */
  2003.     if ( len == 0.0 ) return( 1 ) ;
  2004.     tx /= len ;
  2005.     ty /= len ;
  2006.  
  2007.     /* get dot product for searching */
  2008.     dot = pia->ax * tx + pia->ay * ty ;
  2009.     if ( ( pia->ax * ty > pia->ay * tx ) == pia->flip_edge ) {
  2010.     dot = -2.0 - dot ;
  2011.     }
  2012.  
  2013.     /* binary search through angle list and find matching angle pair */
  2014.     lo = 0 ;
  2015.     hi = pia->hi_start ;
  2016.     while ( lo <= hi ) {
  2017.     ind = (lo+hi)/ 2 ;
  2018.     if ( dot < pia->pis[ind].dot ) {
  2019.         hi = ind - 1 ;
  2020.     } else if ( dot > pia->pis[ind+1].dot ) {
  2021.         lo = ind + 1 ;
  2022.     } else {
  2023.         goto Foundit ;
  2024.     }
  2025.     }
  2026.     /* should never reach here, but just in case... */
  2027.     fprintf( stderr,
  2028.         "Hmmm, something weird happened - bad dot product %lg\n", dot);
  2029.  
  2030.     Foundit:
  2031.  
  2032.     /* test if the point is outside the wedge's exterior edge */
  2033.     pis = &pia->pis[ind] ;
  2034.     inside_flag = ( pis->ex * point[X] + pis->ey * point[Y] <= pis->ec ) ;
  2035.  
  2036.     return( inside_flag ) ;
  2037. }
  2038.  
  2039. void InclusionCleanup( p_inc_anchor )
  2040. pInclusionAnchor p_inc_anchor ;
  2041. {
  2042.     free( p_inc_anchor->pis ) ;
  2043.     free( p_inc_anchor ) ;
  2044. }
  2045. #endif
  2046.